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I.  INTRODUCTION 


The  transient  response  of  a  thin-walled  aluminum  cylinder  subjected  to  an  enveloping 
shock  wave  has  been  simulated  with  the  ADINA  finite  element  program  and  the  results  have 
been  compared  with  counterpart  experimental  data.  The  analysis  employed  the  ADINA  16- 
noded  rectangular  curved  shell  element.  Using  a  linear  analysis,  the  computed  strain  histories 
at  gage  locations  were  found  to  correlate  reasonably  well  with  test  records  for  the  first  4 
milliseconds,  but  some  strain  records  rapidly  diverged  thereafter.  A  subsequent  nonlinear 
analysis  showed  dynamic  buckling  to  be  the  principal  cause  of  this  discrepancy.  The  nonlinear 
calculations  and  test  results  were  found  to  generally  agree  in  the  number  of  buckling  nodes,  the 
magnitude  of  the  deformation,  and  frequencies  excited,  although  buckling  sensitivity  to 
imperfections  prevented  an  exact  coincidence  of  circumferential  modes. 


II.  BACKGROUND 

This  report  presents  a  comparison  between  a  finite  element  analysis  of  the  transient 
response  of  a  thin-walled  aluminum  cylinder  to  an  enveloping  blast  wave  and  the  test  results 
from  a  counterpart  experiment.  The  commercially  available  ADINA1  computer  program  was 
used  for  this  study  because  of  the  variety  of  transient  response  solution  methods  at  its  disposal. 
This  work  is  part  of  an  effort  to  investigate  the  susceptibility  of  lightweight  structures  to  the 
combined  thermal/blast  effects  of  nuclear  explosions  and  to  assess  the  capabilities  of  finite 
element  programs  in  accurately  predicting  the  resulting  structural  response. 

The  basis  for  the  comparison  is  series  of  tests  performed  by  Pearson,  et  al2,  in  1980  in 
which  cylinder  specimens  were  exposed  in  the  BRL  2.2  meter  shock  tube  to  a  blast  only,  thermal 
only,  and  a  combined  thermal/blast  loading.  Thin- walled  cylinders  were  chosen  as 
representative  of  the  lightweight  construction  typical  of  missile  and  airframe  bays.  During  each 
test,  a  cylinder,  with  both  ends  clamped,  was  positioned  with  its  axis  perpendicular  to  the 
direction  of  the  blast  wave  and/or  thermal  radiation.  By  maintaining  the  intensities  of  the 
blast  and  thermal  loadings  constant  for  the  series  of  tests,  a  synergistic  effect  was  conclusively 
demonstrated,  in  that  the  response  to  the  combined  thermal/blast  loading  was  significantly 
greater  than  the  responses  to  the  individual  loadings. 

In  1982,  Gregory  and  Pearson3  attempted  to  reproduce  the  experimental  results  using 
solid  brick  elements  in  the  ADINA  program  to  model  the  cylinder,  but  the  comparison  with 
measured  deflections  was  not  very  good.  At  the  time,  the  poor  correlation  was  attributed  to  an 
excessive  distortional  stiffness  caused  by  the  use  of  single  brick  elements  through  the  shell 
thickness,  although  no  attempt  had  been  made  to  compare  with  the  strain  gage  records. 
Never-the-less,  this  first  attempt  did  suggest  that  perhaps  an  element  more  suitable  to  modeling 
thin-walled  structures  would  perform  better.  For  this  reason,  the  present  analysis  employed 
shell  elements. 

This  investigation  concentrated  on  selecting  an  appropriate  element  to  model  the  thin- 
walled  aluminum  cylinder  of  the  experiment  and  on  evaluating  the  performance  of  the  resulting 
model  against  the  blast  only  results.  One  reason  the  blast  only  test  was  chosen  as  a  basis  of 
comparison  was  that,  of  the  three  loading  combinations,  its  loading  and  response  records  were 
the  cleanest  and  most  complete.  More  specifically,  detailed  transient  blast  loading  data  had 
been  collected  and,  in  addition  to  the  final  deformed  configuration  being  measured,  a  substantial 
number  of  strain  gages  had  been  employed  to  record  the  transient  response.  A  second  reason 
for  the  choice  was  that  the  successful  modeling  of  this  simplest  case  would  provide  a  sound  basis 
for  modeling  the  more  complex  thermal  loading  cases.  Moreover,  a  study  of  the  test  results 
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obtained  by  Pearson,  et  al2,  showed  that  even  this  case  displays  many  of  the  complex 
characteristics,  such  as  buckling,  that  became  more  prominent  in  the  other  two  cases.  In  other 
words,  continuation  of  the  study  would  hinge  on  successfully  modeling  the  blast  only  case. 


m.  CYLINDER  SPECIFICATIONS 

The  test  involved  exposing  a  thin-walled  cylinder  to  a  blast  wave  generated  in  the  BRL 
2.2  meters  diameter  shock  tube4.  The  cylinder  was  positioned  in  the  test  section  of  the  shock 
tube  with  its  axis  perpendicular  to  that  of  the  shock  tube,  so  that  the  blast  wave  enveloped  the 
cylinder  circumferentially.  The  test  specimen  was  manufactured  from  a  rectangular  sheet  of 
6061-T6  aluminum  that  was  rolled  into  a  cylinder  and  butt  welded  along  the  seam.  During  the 
test,  care  was  taken  to  have  this  seam  facing  away  from  the  direction  of  maximum  load 

intensity.  To  achieve  the  clamped  edge  boundary  conditions,  strap  clamps  were  used  to  fix  the 

two  ends  of  the  cylinder  onto  two  relatively  rigid  steel  cylinders,  which,  in  turn  were  attached 
to  the  test  fixture  so  as  to  prevent  the  relative  motion  of  the  two  ends. 

The  specimen  dimensions  were: 

L  —  length  between  clamped  ends  =  0.8  m  (31.5  in.) 

D  —  inside  diameter  =  0.3048  m  (12.0  in.) 

h  =  wall  thickness  =  1.016  mm  (0.04  in.) 

giving  a  D/h  ratio  of  300,  well  within  thin  shell  theory.  The  mechanical  properties  as  obtained 

from  a  tensile  test  were: 

E  —  Young’s  modulus  =  64.73  GPa  (9388  kpi) 

v  —  Poisson’s  ratio  =  0.3285 

yield  stress  =  301  MPa  (43.7  kpi) 

yield  strain  =  0.4195  % 

plastic  tangent  modulus  =  650  MPa  (94.3  kpi) 

The  mass  density  was  taken  as  the  nominal  value  for  aluminum  at  55  degrees  Celsius: 

p  =  density  =  2700  kg/m'  (5.24  slug/ft3) 


IV.  FINITE  ELEMENT  MODEL 

The  finite  element  model  of  the  cylinder  is  depicted  in  Fig.  1.  The  figure  illustrates  the 
global  coordinates  and  the  local  basis  vectors  employed.  By  taking  advantage  of  the  symmetry 
of  the  structure  and  the  assumed  symmetry  of  the  applied  loads,  only  a  quarter  of  the  cylinder 
had  to  be  modeled.  This  was  accomplished  by  placing  symmetry  planes  through  the  cylinder 
axis  and  perpendicular  to  the  axis  midway  from  the  clamped  ends.  Hence,  only  one  clamped 
boundary  was  modeled,  the  upper  boundary  in  Fig.  1. 

The  ADINA  curved  shell  element  with  16  mid-surface  nodes5  was  employed,  with  6 
elements  in  the  axial  direction  and  12  elements  in  the  circumferential  direction,  resulting  in  a  72 
element  modeling  with  703  nodes.  The  choice  of  12  elements  in  the  circumferential  direction 
was  dictated  by  the  location  of  the  pressure  transducers  at  15  degree  intervals  on  the  upper  test 
fixture  cylinder.  In  this  way,  the  element  boundaries  were  made  to  coincide  with  locations  at 
which  pressure  histories  were  experimentally  measured.  In  setting  up  the  problem,  it  was 
necessary  to  not  only  specify  the  nodal  coordinates,  but  also  to  define  the  initial  directions  of 
the  normals  to  the  cylinder,  Vn,  at  the  37  circumferential  nodal  stations,  as  illustrated. 


8 


Figure  1.  ADINA  finite  element  model  of  one  quarter  of  the  cylinder  using  16-noded  shell 
elements. 
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Boundary  conditions  were  specified  by  limiting  the  degrees-of-freedom  (DOF)  along  the 
boundaries:  on  the  z-x  symmetry  plane  no  displacements  in  the  y  direction  and  no  rotations 
about  the  Vt  direction,  on  the  y-z  symmetry  plane  no  displacements  in  the  x  direction  and  no 
rotations  about  the  V2  direction,  and  on  the  clamped  edge  no  displacements  or  rotations  at  all. 
Moreover,  rotation  about  the  normal  vector,  V  ,  were  disregarded  at  all  nodes  because  this  DOF 
is  undefined  for  the  shell  element  employed5.  Consequently,  the  resulting  model  required  3184 
DOF  for  the  703  nodes  employed. 


V.  FREQUENCY  ANALYSIS 

As  a  preliminary  check  on  how  accurately  the  ADINA  model  would  simulate  the 
transient  response  of  the  cylinder,  the  first  20  natural  frequencies  were  calculated  using  the 
curved  shell  element  and  compared  with  predictions  from  a  Rayleigh-Ritz  modal  analysis  of  the 
Donnell  shell  equations  described  in  Kraus6.  The  accuracy  of  the  3-noded  triangular  plate-shell 
was  also  evaluated  by  calculating  the  first  6  frequencies  using  an  equivalent  finite  element  model 
obtained  by  replacing  the  72  curved  shell  elements  by  1296  triangular  plate-shell  elements,  while 
retaining  the  same  nodal  distribution.  Both  analyses  employed  the  subspace  iteration  method  of 
solution,  which  was  found  to  be  more  effective  than  the  determinant  search  method.  Typical  of 
the  results  obtained  is  the  comparison  of  the  first  4  modes  presented  in  Table  1,  where  m  is  the 
number  of  axial  half  waves  and  n  is  the  number  of  circumferential  waves.  In  general,  the 
correlation  between  the  two  finite  element  calculations  and  the  Rayleigh-Ritz  analysis  was 
excellent  over  the  range  of  frequencies  investigated.  From  an  engineering  view  point,  the  results 
using  the  two  finite  element  models  were  essentially  the  same,  but  the  curved  shell  element 
model  was  found  to  required  16  %  fewer  DOF.  For  the  details  of  this  investigation,  including  a 
numerical  convergence  study,  the  interested  reader  may  consult  Santiago,  et  al8. 


Table  1.  Comparison  of  computed  cylinder  frequencies. 


Mode  Shape 

Circular  Frequency  (radian/second) 

m 

n 

Rayleigh-Ritz 

16-Node  Quad 

3-Node  Tri 

i 

4 
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i 

5 
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1894 

1862 

i 

6 

2443 

2464 

2422 

i 

3 

2451 

2395 

2338 

VI.  BLAST  LOADING  MODEL 

The  transient  blast  loading  data  were  taken  from  shot  8-80-7  of  ref.  2.  These 
experimental  data,  in  the  form  of  pressure-time  records,  are  for  a  shock  wave  with  a  peak  static 
overpressure  of  43.7  kPa  (6.34  psi).  As  already  noted,  the  pressure  histories  were  recorded  by  a 
series  of  pressure  gages  mounted  circumferentially  at  15  degree  intervals  on  the  non-deforming 
upper  cylindrical  fixture,  which  formed  a  continuation  of  the  test  cylinder’s  surface.  For  the 
purpose  of  the  present  analysis,  the  pressure  history  at  each  circumferential  location  was 
assumed  to  apply  to  the  entire  cylinder  generator.  In  other  words,  the  spatial  distribution  of 
the  pressure  on  the  cylinder  was  assumed  to  vary  only  circumferentially  and  to  be  uniform  in 
the  axial  direction.  Moreover,  since  the  gages  were  rigidly  mounted,  coupling  between  loading 
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Figure  2.  Pressure-time  loading  functions  derived  from  experimental  data  recorded  at 
circumferential  stations. 
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and  deflection  of  the  specimen  could  not  be  taken  into  account,  but  this  effect  can  be  presumed 
negligible  for  the  relatively  small  deflections  involved. 

The  pressure-time  records  were  digitized  for  the  first  10  milliseconds  as  piece-wise  linear 
histories,  as  illustrated  by  Fig.  2.  Some  of  the  more  prominent  fluid  dynamical  features 
displayed  by  these  graphs  are: 

•  The  increased  delay  in  the  arrival  time  of  the  shock  front  as  the 
angle  $  increases,  with  a  maximum  of  0.95  milliseconds  at  180 
degrees. 

•  The  overpressure  at  the  face-on  location  (zero  degrees)  decaying 
from  105  kPa  normal  reflected  pressure  to  54.4  kPa  stagnation 
pressure. 

•  The  secondary  rise  in  the  pressure  which  emanates  from  the  back 
of  the  cylinder  (detectable  in  the  graphs  from  45  to  150  degrees) 
caused  by  the  shock  fronts  on  each  side  of  the  cylinder  meeting  at 
the  180  degree  location  and  proceeding  to  re-envelop  the  cylinder. 

For  more  details  on  the  shock  envelopment  process,  the  reader  is  directed  to  the  excellent 
pictorial  description  in  Fig.  5  of  ref.  2. 

As  alluded  to  earlier,  the  ADINA  program  accepts  pressure  histories  defined  at  the 
corner  nodes  of  each  element  and  linearly  interpolates  the  corner  values  over  the  element  in 
order  to  calculate  the  equivalent  nodal  forces.  The  interpolation  is  accomplished  by  defining 
over  each  set  of  four  elements  sharing  a  common  node  a  pyramid  function  that  assumes  the 
value  of  the  node  pressure  at  the  apex  node  and  that  vanishes  at  the  corner  nodes.  The 
pressure  over  an  element  is  then  easily  found  by  summing  the  contributions  from  the  pyramid 
functions  at  the  four  corner  nodes.  This  procedure  is  illustrated  in  Fig.  3  for  the  present  (two- 
dimensional)  case,  in  which  dash  lines  represent  the  pyramid  functions  defined  for  each  corner 
node  and  the  solid  line  represents  the  full  pressure  acting  on  each  element,  obtained  by 
summing  the  contributions  from  the  corner  node  pyramid  functions.  The  program  then  uses  the 
resultant  pressure  to  compute  the  equivalent  forces  at  all  nodes,  both  corner  and  interior,  in  a 
consistent  manner. 


VH.  STRAIN  CALCULATIONS 

In  performing  an  analysis,  the  ADINA  program  computes  strains  at  the  Gauss 
integration  points  of  each  element.  However,  only  the  values  at  the  corner  Gauss  points  are 
easily  retrievable  for  comparison  with  experimental  data.  The  present  analysis  used  the  default 
values  of  3X3X2  for  the  Gauss  point  distribution  associated  with  each  16-noded  shell  element, 
resulting  in  the  distribution  depicted  in  Fig.  4,  wherein  cell  coordinates  are  used  to  specify  point 
locations  and  the  open  circles  denote  the  points  at  which  strain  values  are  saved. 

During  the  experiment,  a  total  of  ten  strain  gages  placed  on  the  inner  surface  of  the 
cylinder  were  used  to  record  transient  strain  histories.  Fig.  5  shows  the  locations  and  directions 
of  these  gages  superimposed  on  a  developed  projection  of  the  finite  element  model.  This  figure 
assumes  the  symmetric  deformation  of  the  cylinder,  and,  hence  does  not  distinguish  between 
gage  locations  on  the  upper  or  lower  and  left  or  right  quarters  of  the  cylinder  The  precise 
locations  of  these  gages  in  terms  of  the  axial  distance  from  the  mid-plane  and  circumferential 
angle  are  presented  in  Table  2. 
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Figure  3.  Generation  of  the  pressure  over  the  interior  of  an  element  by  a  linear  interpolation 
of  the  pressure  functions  defined  at  corner  nodes. 
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Figure  4.  Location  of  the  Gauss  integration  points  employed  in  a  typical  16-noded  shell 
element. 
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Figure  5. 


Geometric  locations  and  orientations  of  the  strain  gages  on  the  cylinder  and  the 
locations  of  the  adjacent  Gauss  points  used  for  interpolation. 
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Table  2.  Locations  of  strain  gages  on  test  cylinder. 


Gage 

number 

Gage 

direction 

Angle  from  face-on 
position  (degree) 

Distance  from 
mid-plane  (cm) 

1 

axial  &  circ. 

0 

0.0 

2 

axial  &  circ. 

0 

19.7 

3 

axial 

0 

-37.5 

4 

circumferential 

315 

0.0 

5 

axial  &  circ. 

270 

0.0 

6 

axial 

270 

-37.5 

7 

circumferential 

135 

0.0 

Since  the  gages  recorded  strains  at  the  inner  surface,  while  the  program  computed 
strains  at  Gauss  points  located  in  the  shell  interior  (see  Fig.  4),  it  was  necessary  to  linearly 
extrapolate  interior  values  along  the  normal  or  t  direction  to  the  inner  surface.  Moreover,  as 
Fig.  5  makes  clear,  the  geometric  locations  of  the  gages  do  not  coincide  with  the  Gauss 
integration  points  at  which  the  computed  strains  were  stored,  so  that  some  manipulation  of  the 
computed  values  was  required  to  compare  with  experimental  records.  This  was  accomplished  by 
a  simple  linear  interpolation  of  the  strain  values  computed  at  the  adjacent  Gauss  points;  Fig.  5 
shows  the  Gauss  points  adjacent  to  each  gage  location  as  open  circles.  Notice  that  only  two 
Gauss  points  were  required  for  gage  locations  on  the  symmetry  boundaries  and  only  one  for 
those  at  the  intersection  of  symmetry  boundaries,  since  there  are  mirror  image  Gauss  points 
symmetrically  located  outside  the  boundaries. 

Because  the  program  calculates  strain  components  relative  to  global  xyz  coordinates,  see 
Fig.  1,  rather  than  shell  tangent  coordinates,  a  transformation  was  necessary  in  order  to  obtain 
components  in  the  gage  directions  at  some  locations.  For  the  axial  strains  at  gage  locations  1,  2, 
3,  5  and  6  in  Fig.  5,  no  transformations  were  required,  since  they  are  given  directly  by  the  eyy 
components.  Neither  were  transformations  of  the  circumferential  strains  at  gage  locations  1  and 
2  necessary,  since  they  are  considered  sufficiently  well  approximated  by  exx1  nor  for  that  at  gage 
location  5  since  ezz  is  sufficiently  close  there.  Only  at  gage  locations  4  and  7  were 
transformations  necessary  in  order  to  determine  the  circumferential  strains.  The  transformation 
was  accomplished  easily  by  applying  Mohr’s  circle  formula9: 

ecir  =  7  lezz  +  -{ezz-  ezx)  cos  29  -  exz  sin  29}  (1) 


where  6  is  the  angle  between  the  z-axis  and  the  normal  vector,  see  Fig.  1,  The  circumferential 
strains  at  locations  4  and  7  were  then  obtained  by  simply  substituting  45  and  135  degrees  for  6 
in  this  formula. 

In  actuality,  it  was  found  more  effective  to  program  the  above  procedure  in  somewhat 
reverse  order: 

•  First,  where  necessary,  (1)  transformed  the  strain  components  at 
the  Gauss  points  adjacent  to  gage  locations  of  interest  to 
circumferential  strains. 

•  Second,  strain  values  at  all  Gauss  points  neighboring  gage 
locations  were  extrapolated  to  the  inner  surface. 
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•  Third,  the  just  computed  surface  values  adjacent  to  each  gage 
location  were  interpolated  to  obtain  the  computed  values  of  the 
strains  in  the  appropriate  directions  (axial  or  circumferential)  at 
the  gage  locations. 

In  addition  to  permitting  the  computed  strains  to  be  compared  with  the  experimental  strains, 
this  procedure  provided  some  degree  of  smoothing  of  the  strains  across  element  boundaries. 


VIIL  LINEAR  ANALYSIS 

A  linear  analysis  was  first  performed  to  determine  the  strain  history  at  each  gage 
location  for  comparison  with  the  experimental  records.  The  elastic  material  properties, 
summarized  in  Chapter  III,  were  read  in  as  initial  input  data  The  loading  was  supplied  at  each 
time  step  from  the  pressure  histories  as  described  in  Chapter  VI.  As  in  the  frequency  analysis,  a 
consistent  mass  matrix  formulation  was  employed.  The  implicit  (Newmark)  time  integration 
method  was  used  with  a  —  0.25  and  b  —  0.5,  (the  so-called  trapezoid  rule).  The  calculation 
proceeded  for  200  time  steps  using  an  increment  of  50  microseconds,  for  a  total  real  time  of  10 
milliseconds. 

The  increment  of  50  microseconds,  chosen  for  convenience,  was  sufficiently  small  to 
adequately  model  the  first  20  vibratory  modes  predicted  by  the  frequency  analysis,  since  the 
highest  mode  calculated  had  a  period  of  847  microseconds.  Moreover,  the  test  records,  Fig.  6 
and  7,  indicate  that  the  period  of  the  dominant  response  frequency  was  in  the  order  of  1.7 
milliseconds,  more  than  an  order  of  magnitude  greater  than  the  time  increment. 

The  comparison  between  computed  strain  histories  and  the  experimental  records  is 
shown  in  Fig.  6  and  7.  The  data  compare  rather  well  for  the  axially  oriented  gages  (Fig.  6).  For 
instance,  the  comparisons  of  frequencies  at  locations  5  and  6,  as  well  as  the  first  4  to  5  ms  at 
location  1,  show  excellent  agreement.  Only  a  sight  phase  shift  is  apparent,  which  can  be 
attributed  to  stiffening  due  to  the  use  of  a  consistent  mass  formulation,  see  Strang  and  Fix10 
Moreover,  except  at  location  3  and  the  last  5  milliseconds  at  location  1,  relative  differences  in 
the  amplitudes  are  very  small. 

It  is  interesting  to  contrast  the  correlations  obtained  at  the  two  gage  locations  near  the 
boundary;  while  that  at  location  6  is  excellent,  the  one  at  location  3  is  not  as  good  This  implies 
that  the  clamped  edge  condition  was  satisfied  better  at  the  location  side-on  to  the  blast  wave 
than  at  the  face-on  location,  where  some  slipping  may  have  occurred.  This  is  not  too  surprising, 
since  the  face-on  location  received  the  full  brunt  of  the  reflected  pressure  and,  hence,  the  edge  is 
more  apt  to  slip  there. 

The  comparisons  between  the  computed  and  recorded  strain  histories  in  the 
circumferential  direction  are  noticeably  poorer  (Fig.  7).  Even  during  initial  times,  frequencies  do 
not  correlate  nearly  as  well.  Moreover,  with  the  exception  of  location  4,  the  experimental  records 
begin  to  diverge  noticeably  from  the  computed  results  after  3  to  5  milliseconds  Though  not  as 
pronounced,  this  unusual  behavior  can  also  be  detected  in  the  axial  histories  at  locations  I  and 
possibly  3,  although  the  erratic  correlation  at  the  latter  location  makes  a  definitive  statement 
difficult.  The  records  of  the  circumferential  strains  along  the  mid-section  of  the  cylinder 
(locations  1,  4,  5,  and  7)  imply  that  the  divergent  behavior  may  be  associated  with  a  pattern  of 
circumferential  waves,  suggesting,  perhaps,  that  some  buckling  has  occurred.  This  conjecture  is 
confirmed  somewhat  by  the  noticeably  smaller  effect  on  the  axial  records.  To  explore  the 
possibility  of  buckling  using  ADINA,  a  nonlinear  analysis  was  necessary. 
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Comparison  of  linearly  computed  and  experimental  strain  histories  at  axially 
oriented  gage  locations. 
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Figure  7.  Comparison  of  linearly  computed  and  experimental  strain  histories  at 
circumferentially  oriented  gage  locations. 
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IX.  BUCKLING  LOAD  CALCULATION 


As  a  first  step  the  critical  buckling  pressure  and  corresponding  mode  shape  were 
determined  using  the  von  Mises  equation11  for  a  cylinder  subjected  to  a  uniform  external 
pressure: 


-  M  ( 4!  +  _L_(  A)  V2-  1  +  (2„2  - 1  -  ,)4 

a  l  n2-l  1 2(1— i^2)  l  aj  I 


K  = 


1  +  ( nL/xa, )2 


(2) 

(3) 


where  a  (=  .1529  m)  is  the  mean  radius  of  the  cylinder,  n  is  the  number  of  circumferential 
waves,  and  the  other  parameters  have  been  previously  defined  in  Chapter  III.  Although  this 
equation  applies  strictly  to  hinged  boundary  conditions  only,  because  the  length  to  radius  ratio 
is  greater  than  5,  the  boundaries  can  be  considered  remote  enough  to  allow  its  application  to 
clamped  edge  boundaries.  Substituting  the  values  for  the  test  cylinder  in  this  expression  yields 
a  threshold  pressure  of  42  kPa  for  an  associated  buckling  mode  with  4  circumferential  waves. 

Although  the  experimentally  measured  pressures,  as  shown  in  Fig.  2,  are  not  spatially 
uniform,  we  do  observe  that  the  pressures  at  the  first  three  stations  decay  to  approximately  60 
kPa  after  their  initial  peaks  of  approximately  lOOkPa,  while  those  at  the  remaining  stations 
hover  at  about  40  kPa.  Hence,  the  experimental  pressures  can  be  expected  to  induce  buckling, 
the  average  clearly  being  in  the  range  of  the  predicted  critical  buckling  pressure.  Moreover,  if 
the  buckling  did  indeed  occur  dynamically,  the  theory  of  dynamic  buckling  would  predict  that 
mode  shapes  corresponding  to  the  locally  higher  pressures  are  more  likely  to  be  observed12. 
Solving  (2)  and  (3)  for  the  next  higher  modes,  we  find  a  pressure  of  47  kPa  corresponding  to  5 
circumferential  wave  and  65  kPa  corresponding  to  6  circumferential  waves,  well  within  the 
range  of  the  experimentally  observed  pressures  and  implying  that  5  and  possibly  even  6 
circumferential  waves  can  be  expected  to  occur.  The  post  test  measurement  of  the  cylinder 
cross-section,  Fig.  8,  appears  to  verify  that  possibility,  implying  that  perhaps  five 
circumferential  waves  were  involved. 


X.  NONLINEAR  ANALYSIS 

With  this  evidence  that  the  cylinder  may  have  experienced  some  buckling,  a  geometric 
nonlinear  transient  analysis  was  undertaken  with  the  same  loading  functions  as  used  in  the 
linear  analysis.  The  nonlinear  analysis  employed  the  total  Lagrangian  formulation  available 
with  the  shell  element.  As  in  the  linear  analysis,  the  trapezoid  rule  was  used  for  implicit  time 
integration,  with  the  same  time  step  of  50  microseconds.  However,  it  was  necessary  to  use  a 
lumped  mass  matrix,  rather  than  a  consistent  mass  matrix,  because  only  the  former  option  has 
been  programmed  into  the  nonlinear  solution  algorithm  for  the  curved  shell  element.  Hence,  the 
nonlinear  analysis  was  not  the  exact  analog  of  the  linear  case. 

The  results  of  the  nonlinear  analysis  confirmed  to  a  great  extent  that  buckling  had  taken 
place.  Indeed,  the  circumferential  strain  histories  computed  at  15  degree  intervals  along  the 
inner  surface  at  the  mid-cylinder  circumference,  as  shown  in  Fig.  9,  reveal  a  significant  increase 
in  their  amplitudes  beginning  at  approximately  4  milliseconds,  corresponding  to  the  onset  of 
buckling  experimentally  observed  in  the  strain  gage  records.  The  computed  mid-cylinder  radial 
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Figure  8.  Post  test  measurements  of  the  mid-cylinder  radial  deflection. 
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Figure  9. 


Circumferential  strain  histories  computed  from  the  nonlinear  solution  on  the  inner 
surface  at  mid-cylinder  stations  15  degree  apart. 
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deflections,  depicted  in  Fig.  10,  corroborate  the  last  result  by  also  exhibiting  a  noticeable 
increase  after  4  milliseconds.  Moreover,  the  deflected  cross-sections,  despite  some  erratic 
oscillations,  indicate  a  fairly  persistent  pattern  of  five  circumferential  waves,  with  the  locations 
of  the  maxima  and  minima  remaining  fairly  stationary,  in  agreement  with  dynamic  buckling 
theory  and  the  experimental  post-shot  shape  (Fig.  8).  Isometric  computer  plots  of  the  deformed 
surface  of  the  cylinder,  collected  at  select  times  in  the  appendix  at  the  end  of  this  report,  also 
verify  the  occurrence  of  dynamic  buckling.  As  is  to  be  expected,  they  show  that  the 
deformation  is  considerably  more  complex  than  predicted  by  the  buckling  analysis.  For 
example,  from  these  plots  we  see  how  the  longitudinal  deformation  progessively  becomes  more 
complex  with  time:  starting  as  a  single  wave  and  by  6  millisecond  beginning  to  show  the 
development  of  three  waves. 

When  the  computed  strain  histories  are  directly  compared  with  the  corresponding  test 
data  for  the  strain  gages  along  the  mid-cylinder,  Fig.  11  and  12,  the  correlations  are  not  entirely 
satisfactory.  The  computed  histories  for  the  axially  oriented  gages  (Fig  11),  though  somewhat 
choppier  than  those  from  the  linear  solution  (Fig.  6),  are  in  general  agreement  with  these 
computations  and  compare  reasonably  well  with  test  results.  The  computed  histories  for  the 
circumferentially  oriented  gages  (Fig.  12),  although  an  improvement  over  the  linear  solutions 
(Fig.  7)  in  clearly  simulating  the  increase  in  amplitudes  following  buckling,  do  not  correlate  that 
well  with  test  results  at  all  gage  locations.  For  example,  the  computed  strains  at  gage  locations 
1,  5,  and  7  are  certainly  a  considerable  improvement  over  the  linearly  computed  strains  in 
modeling  the  magnitude  and  frequency  during  buckling,  but  the  general  trends  are  not  exactly 
reproduced.  At  gage  location  4  the  linear  analysis  is  clearly  closer  to  the  test  result,  since  the 
gage  did  not  record  an  increased  strain  accompanying  buckling. 

On  the  other  hand,  considering  the  lack  of  exactitude  that  occurs  in  buckling 
phenomena  due  to  uncertainties  in  loading  and  material  properties,  we  thought  it  worthwhile  to 
attempt  to  compare  the  test  records  with  the  computed  strains  at  other  circumferential 
locations  in  order  to  improve  the  correlation,  as  depicted  in  Fig.  13.  Except  at  gage  5,  where 
the  best  correlation  is  still  obtained  with  the  strains  computed  at  90  degrees,  we  notice  an 
immediate  improvement  in  the  correlations.  For  example,  the  correlations  between  the  gage  1 
record  and  the  computed  history  at  135  degree,  and  the  gage  4  record  and  the  15  degrees 
history  show  considerable  improvement  over  those  in  Fig.  12.  And,  indeed,  the  correlation 
between  the  gage  7  record  and  the  150  degrees  history  is  extraordinarily  good.  Ideally,  it  would 
have  been  more  convincing  if  a  consistent  (and  hopefully  small)  rotation  of  the  computed  strain 
stations  about  the  cylinder  axis  had  resulted  in  an  improved  correlation  at  all  the  strain  gage 
locations,  but,  despite  a  number  of  attempts,  no  such  rotation  was  discovered. 

It  is  important  to  point  out  that  the  computed  strains,  in  agreement  with  experimentally 
observations  (Fig.  7),  never  exceeded  the  plastic  yield  strain,  implying  that  at  least  for  the  first 
10  milliseconds  only  elastic  buckling  has  taken  place.  To  make  sure  that  the  use  of  a  purely 
elastic  analysis  was  valid,  an  additional  geometric  and  material  nonlinear  analysis  was 
performed  and  was  found  to  give  identical  results.  The  elastic  behavior  of  the  deformation  to 
some  extent  explains  the  oscillation  observable  in  Fig.  9,  since  there  was  no  plastic  dissipation  to 
absorb  the  kinetic  energy  and,  hence,  “freeze”  the  buckled  pattern. 

However,  the  conclusion  that  only  elastic  deformation  took  place,  does  somewhat  conflict 
with  the  post-test  measurement  of  permanent  deformation  in  the  cross-section,  Fig.  8. 
Accepting  that  plastic  yielding  most  likely  occurred  after  10  milliseconds,  there  is  the  question 
of  what  caused  this  to  happen.  Possible  causes  are: 

•  Intial  imperfections  in  the  geometry  of  the  test  specimen  (a  certain 
amount  of  out-of-roundness  was  detected  in  specimens2). 
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Figure  10. 


Deflection  pattern  of  the  mid-cylinder  cross-section  at  select  times  (in  milliseconds) 
as  computed  by  the  nonlinear  analysis. 
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Figure  11. 


Comparison  of  axial  strain  histories  recorded  at  mid-cylinder  gage  locations  and 
computed  by  the  nonlinear  analysis  at  corresponding  stations. 
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Figure  12.  Comparison  of  circumferential  strain  histories  recorded  at  mid-cylinder  gage 
locations  and  computed  by  the  nonlinear  analysis  at  corresponding  stations. 
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Figure  13.  Comparison  of  circumferential  strain  histories  recorded  at  mid-cylinder  gage 
locations  and  computed  by  the  nonlinear  analysis  at  circumferential  stations  selected 
for  optimum  correlation. 
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•  Uncertainties  and  lack  of  uniformity  in  the  material  properties  of 
the  experimental  cylinder,  especially  the  plastic  yield  stress  (recall 
the  welded  seam  and  some  yielding  caused  by  using  strap  clamps 
at  the  ends). 

•  The  vibratory  modes  excited  by  the  nonuniform  loading 
eventually  combining  to  produce  a  deformation  that  exceeded  the 
yield  strain  at  some  points. 

•  Uncertainty  about  the  pressure  loading  distribution  during  the 
first  10  milliseconds  (recall  the  assumption  of  axial  uniformity). 

•  Changes  in  the  intensity  and  distribution  of  the  pressure  loading 
after  the  first  10  milliseconds. 

Of  these  possibilities,  it  is  most  unlikely  that  geometric  imperfections  or  local  weaknesses  in  the 
material  caused  plastic  deformation,  since  their  effects  would  have  been  manifested  soon  after 
buckling  occurred  during  the  first  10  milliseconds.  More  likely,  but  not  too  probable  for  the 
same  reason,  is  the  possibility  that  the  vibratory  modes  eventually  combined  to  produce 
plasticity.  The  assumption  of  axial  uniformity  of  the  pressure  distribution  is  probably  quite 
good  for  the  first  10  milliseconds,  giving  a  fairly  close  simulation  of  loading  imposed  on  the 
cylinder,  since  not  enough  time  had  elapsed  for  reflections  from  the  test  fixture  and  the  ends  of 
the  shock  tube  to  distort  the  flow  in  the  test  section.  Hence,  the  likeliest  reason  for  the 
permanent  plastic  deformation  is  the  changes  that  the  pressure  loading  undergoes  during  late- 
times.  More  specifically,  as  the  air  emptied  out  of  the  shock  tube,  rarefaction  waves  proceeding 
from  the  open  end  probably  accelerated  the  flow  velocity  sufficiently  to  distort  the  distribution 
of  the  drag  loads  on  the  cylinder  and  locally  magnify  the  pressure  level,  resulting  in  the  further 
(plastic)  deformation  of  the  cylinder.  Moreover,  the  resulting  late-time  pressure  distribution 
would  probably  have  been  axially  nonuniform,  making  the  gage  readings,  had  they  been 
recorded  at  all,  unrepresentative  of  the  loading  on  the  cylinder.  For  this  reason,  the  lack  of 
late-time  loading  data,  and  the  extensive  computational  effort  that  it  would  have  involved,  no 
attempt  was  made  to  carry  the  calculations  past  the  initial  10  milliseconds. 

Whatever  the  cause  of  the  plastic  deformation,  it  is  significant  that  both  the  von  Mises 
buckling  analysis  and  the  ADINA  nonlinear  solution  predict  a  buckling  pattern  very  similar  to 
that  measured  on  the  cylinder.  Moreover,  since  it  is  likely  that  plastic  yielding  would  occur 
where  the  strains  are  extreme  (i.e.,  where  the  curvature  is  greatest)  and  since  the  computed 
buckling  pattern,  although  elastic,  does  not  change  much  with  time  (see  Fig.  10),  it  is  reasonable 
to  expect  that  the  same  pattern  would  persist  into  the  late-time  plastic  flow  phase,  justifying  to 
a  certain  degree  the  comparison  of  the  post  test  shape  and  the  computed  elastic  cross-sections. 


XI.  SUMMARY  AND  CONCLUSIONS 

The  subspace  iteration  technique  in  the  ADINA  code  was  found  to  effectively  calculate 
the  frequencies  and  mode  shapes  of  the  test  cylinder.  The  first  20  modes  were  calculated  using 
the  ADINA  16-noded  shell  element  and  the  first  6  modes  using  the  ADINA  3-noded  triangular 
plate  element.  Frequencies  and  mode  shapes  predicted  by  both  elements  were  found  to  correlate 
excellently  with  each  other  and  with  those  from  a  Rayleigh-Ritz  analysis  based  on  the  Donnell 
shell  equations. 

Good  overall  correlation  with  the  data  from  the  shock  tube  cylinder  test  was  obtained 
using  the  ADINA  16-noded  curved  shell  element.  The  strain  histories  calculated  from  the  linear 
solution  corresponded  fairly  well  with  the  strains  recorded  for  the  axially  oriented  gages  over  the 
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entire  interval.  The  correspondence  with  the  strains  for  the  circumferentially  oriented  gages 
was  also  good,  but  only  for  the  first  4  milliseconds,  after  which  the  strains  at  some  gage 
locations  show  an  increasing  divergence.  The  subsequent  geometric  nonlinear  analysis  showed 
dynamic  buckling  to  be  the  principal  cause  of  this  divergence.  The  strain  histories  computed  by 
the  nonlinear  analysis  showed  the  same  significant  increase  in  the  circumferential  components 
occurring  after  4  milliseconds.  Although  the  circumferential  strain  gage  data  did  not  completely 
agree  with  the  computed  stain  histories  at  the  corresponding  gage  locations,  other  locations  were 
found  at  which  agreement  was  considerably  improved. 

The  deflection  pattern  predicted  by  the  nonlinear  analysis  and  the  post-test  deflection 
measurements  were  found  to  agree  substantially  in  the  number  of  buckling  waves  and  the 
magnitude  of  the  deformation.  That  the  buckling  pattern  on  the  cylinder  did  not  coincide 
exactly  with  that  measured  on  the  test  specimen  was  most  probably  due  to  geometric  and 
material  imperfections  and  the  sensitivity  of  buckling  phenomena  to  such  defects.  However, 
given  the  difficulty  of  making  accurate  buckling  predictions  due  to  uncertainties  in  material 
properties,  specimen  geometry,  and  loads  definition,  and  that  perhaps  a  less  crude  finite  element 
mesh  should  have  been  used,  the  finite  element  analysis  succeeded  in  generally  reproducing  the 
salient  features  of  measured  deformation  to  a  remarkable  degree. 

This  investigation  conclusively  demonstrated  that  dynamic  buckling  had  a  significant 
effect  on  the  blast  response  of  the  cylinder  and  that  the  ability  of  the  finite  element  calculations 
to  accurately  predict  the  deflections  hinged  on  correctly  modeling  this  phenomenon.  The 
calculations  also  emphasized  the  importance  of  choosing  an  appropriate  finite  element  modeling 
to  perform  the  transient  analysis  of  shell-like  structures,  in  particular,  whenever  there  is  a 
reasonable  likelihood  of  dynamic  buckling  occurring. 
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APPENDIX:  CYLINDER  SURFACE  DEFLECTIONS 


The  following  pages  of  this  appendix  contain  computer  plots  of  the  deformed  surface  of 
the  cylinder  at  select  times.  They  were  calculated  by  the  ADINA  program  using  the  total 
Lagrangian  method  for  modeling  geometric  nonlinear  effects.  The  material  was  assumed  to  be 
linearly  elastic.  For  the  real  times  indicated,  each  page  shows  a  longitudinal  cross-section 
through  the  crown  line,  a  lateral  cross-section  through  the  mid-cylinder,  and  an  isometric  of  the 
surface  deflection.  The  deflections  are  shown  magnified  by  a  factor  of  ten  (10)  relative  to  the 
cylinder  dimensions.  The  plots  trace  the  development  of  dynamic  buckling  beginning  at  3.5 
millisecond.  The  deflections  do  not  become  larger  because  no  plastic  yielding  occurs  and  hence 
no  permanent  deformation  can  result. 
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